//
// Configuration.cxx
//
// Created by TB on 11/17/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <cmath>
#include <cstdlib>
#include <sstream>
#include <stdexcept>

#include <getopt.h>

#include "tppmc/Configuration.h"
#include "tppmc/Enumerations.h"

namespace tppmc {
   
   const double PYTHIA_MIN_PARTON_PT(2.5);
   
   enum MachineType{RHIC_200,RHIC_500,FNAL_E704_PROTON,FNAL_E704_PBAR};
   enum GluonSiversType{NO_SIVERS,AS_UP_QUARK,POS_LIMIT,AS_UP_QUARK_POS_BND};
   enum PythiaEventType{QCD2x2,GENJPSI};
   
   Configuration::Configuration(int argc, char* argv[])
   : mPytune(0)
   , mLhapdf(-1, -1) {
      fftype  = 4;  
      order   = LO; 
      numEvts = 0; 
      machine_type = RHIC_200; 
      pyth_evt_type = QCD2x2;
      quark_transv_soffer = 0; 
      quark_siv_posbnd = 0; 
      etamin = -999.;
      etamax = 999.;
      mean_jt2 = 0.25;
      min_hadron_pt = 0.;
      fav_pos_bnd = 0; 
      unfav_pos_bnd = 0; 
      helpflag = 0; 
      sivers_type = 0;
      pyth_parton_pt = PYTHIA_MIN_PARTON_PT;
      primordial_kt = std::sqrt(0.25);
      // Default process selection for QCD 2 --> 2 scattering
      int processes[7] = {11, 12, 13, 28, 53, 68, 96};
      mProcesses = std::set<int>(processes, processes + 7);
      if(argv) {
         Set(argc, argv);
      } // if
   }

   Configuration::~Configuration() {
   }

   int Configuration::Set(int argc, char** argv) {
      /* no_argument, required_argument, optind, option are all defined
       in getopt.h.
       struct option {
          char* name;
          int has_arg;
          int* flag;
          int val;
       };
       */
      option long_options[] = {
         {"help",                   no_argument, &helpflag, 1 },
         {"quark_transv_soffer",    no_argument, &quark_transv_soffer, 1 },
         {"tppmc_frag_function",    required_argument, NULL, 'g'},
         {"ffunc",                  required_argument, NULL, 'f'},
         {"output",                 required_argument, NULL, 'o'},
         {"order",                  required_argument, NULL, 'r'},
         {"events",                 required_argument, NULL, 'e'},
         {"rhic200",                no_argument, &machine_type, RHIC_200 },
         {"rhic500",                no_argument, &machine_type, RHIC_500 },
         {"fnalftp",                no_argument, &machine_type, FNAL_E704_PROTON },
         {"fnalftpb",               no_argument, &machine_type, FNAL_E704_PBAR },
         {"qcd2x2",                 no_argument, &pyth_evt_type, QCD2x2},
         {"processes",              required_argument, NULL, 'P'},
         {"min_hard_pt",            required_argument, NULL, 'm'},
         {"jpsi",                   no_argument, &pyth_evt_type, GENJPSI},
         {"pytune",                 required_argument, NULL, 'p'},
         {"qa",                     required_argument, NULL, 'q'},
         {"lhapdf",                 required_argument, NULL, 'l'},
         {"sivers",                 required_argument, NULL, 's'},
         {"eta",                    required_argument, NULL, 'E'},
         {"meanjt2",                required_argument, NULL, 'j'},
         {"min_hadron_pt",          required_argument, NULL, 't'},
         {"primordial_kt",          required_argument, NULL, 'k'},
         {NULL, 0, NULL, 0}
      };
      while(true){
         int option_index(0);
         int c = getopt_long(argc, argv,"g:f:o:r:e:p:P:l:q:s:E:j:t:m:k:",
                             long_options, &option_index);
         if(c == -1) {
            break;
         } // if
         std::stringstream ss; // For string --> number conversions
         // Process short options
         switch(c) {
            case 0:
               if(long_options[option_index].flag!=0) break;
               printf ("option %s", long_options[option_index].name);
               if (optarg) printf (" with arg %s", optarg);
               printf ("\n");
               break;
            case 'f':
               fftype = strtol(optarg,NULL,10);
               break;
            case 'o':
               outfile = optarg; 
               break;
            case 'r':
               order = strtol(optarg,NULL,10);
               break;
            case 'e':
               numEvts = strtol(optarg,NULL,10);
               break;
            case 'k':
               primordial_kt = strtof(optarg, NULL);
               break;
            case 'p':
               ss.str("");
               ss.clear();
               std::cout << "User passed pythia tune " << optarg << std::endl;
               ss << optarg;
               ss >> mPytune;
               break;
            case 'P':
               ss.str("");
               ss.clear();
               ss << optarg;
               mProcesses.clear();
               while(ss.good()) {
                  int i;
                  ss >> i;
                  mProcesses.insert(i);
               } // while
               break;
            case 'l':
#ifdef USING_LHAPDF
               // Input format is "set,member"
               std::cout << optarg << std::endl;
               ss.str("");
               ss.clear();
               ss << optarg;
               ss >> mLhapdf.first >> mLhapdf.second;
#else
               std::cerr <<
               "LHAPDF functionality is not supported." <<
               "\nRun configure --help for how to compile with LHAPDF support"
               << std::endl;
#endif
               break;
            case 'q':
               ss.str("");
               ss.clear();
               ss << optarg;
               ss >> mQaFiguresFileName;
               break;
            case 's':
               sivers_type = strtol(optarg, NULL, 10);
               break;
            case 'E':
               ss.str("");
               ss.clear();
               ss << optarg;
               ss >> etamin >> etamax;
               break;
            case 'j':
               ss.str("");
               ss.clear();
               ss << optarg;
               ss >> mean_jt2;
               break;
            case 't':
               std::cout << optarg << std::endl;
               ss.str("");
               ss.clear();
               ss << optarg;
               ss >> min_hadron_pt;
               std::cout << min_hadron_pt << std::endl;
               break;
            case 'm':
               pyth_parton_pt = strtof(optarg,NULL); 
               break; 
            default:
               abort();
         } // switch
      } // while
      
      if(optind < argc or 1 == helpflag) {
         PrintUsage(std::cout);
      } // if
      
      if((0 == numEvts or outfile.empty()) and mQaFiguresFileName.empty()) {
         std::cout <<
         "Must specify at least an output file (-o) and number of events (-e) to run." << std::endl;
         exit(-1); 
      } // if
      
      if(not Verify())
         throw std::runtime_error("Invalid Configuration");
      
      return optind;
   }
   
   
   void Configuration::PrintUsage(std::ostream& o) {
      o << std::endl; 
      o << "usage: tppmc -e [numevents] -o [output file]\n";
      o << std::endl; 
      o << "Distribution functions:\n";
      o << "--lhapdf, -l=\"set member\" : LHAPDF set and member " <<
           "for the unpolarised PDF (default: GRV98)\n";
      o << "--sivers=[n] : Sivers distribution function\n" <<
           "     (default) 0 = None\n" <<
           "               1 = Torino group 2009 fit\n" <<
           "               2 = positivity bound from unpolaried PDF\n";
      o << "--quark_transv_soffer : u/d quark transversity at Soffer bound\n";
      o << "--ffunc, -f=[n] : quark and gluon fragmentation type\n";
      o << "                  1 = DSS\n";
      o << "                  2 = Fragment parton into a single hadron\n";
      o << "        (default) 4 = DSS + Collins (pion) Torino group fit\n";
      o << "                  5 = DSS + Collins (pion), positivity bound\n";
      o << "--order, -o=[n] : order\n";
      o << "                  0 = no initial/final-state radiation, " <<
                             "leading order fragmentation functions\n" <<
           "        (default) 1 = initial/final-state radiation, " <<
                             " NLO fragmentation functions\n";
      o << std::endl; 
      o << "Physics parameters:\n";
      o << "--min_hard_pt, -m : minimum pT in pythia hard scattering " <<
           "(default = " << PYTHIA_MIN_PARTON_PT << " GeV/c)\n";
      o << "--eta, -E=\"min max\" : set pseudorapidity range for hadrons\n";
      o << "--min_hadron_pt=[pt], -t : minimum hadron pT (default 0)\n";
      o << "--primordial_kt=[kt], -k : primordial kT of partons in hadrons " <<
           "(default 0.5 GeV/c)\n";
      o << std::endl; 
      o << "Event generation (pick one):\n";
      o << "--qcd2x2 : generate QCD 2x2 scattering (default)\n";
      o << "--jpsi : J/Psi production\n";
      o << "--pytune, -p=[n] : select PYTHIA tune\n";
      o << std::endl; 
      o << "Output:\n";
      o << "--events, -e : number of scattering events to generate\n";
      o << "--output, -o=[filename] : root histogram and output file\n";
      o << "--qa, -q=[file name] : generate QA figures in a ROOT file\n";
      o << "--help : print this help screen\n";
      o << std::endl; 
      o << "Machine Parameters:\n";
      o << "--rhic200 : RHIC collider 200 GeV (default)\n";
      o << "--rhic500 : RHIC collider 500 GeV\n";
      o << "--fnalftp : FNAL E704 p+p fixed target at 200 GeV " <<
           "(uses equivalent CMS energy)\n";
      o << "--fnalftpb : FNAL E704 pbar+p fixed target at 200 GeV " <<
           "(uses equivalent CMS energy)\n";
      o << std::endl; 
   }

   void Configuration::Speak(std::ostream& o) const {
      o << "Output file: " << outfile << std::endl;
      o << "Fragmentation function type: " << fftype << std::endl;
      o << "Fragmentation order: " << order << std::endl;
      o << "Number of events: " << numEvts << std::endl;
      o << "Machine type: " << machine_type << std::endl;
      o << "PYTHIA event type: " << pyth_evt_type << std::endl;
      o << "Sivers set: " << sivers_type << std::endl;
      o << "Minimum hard pT: " << pyth_parton_pt << std::endl;
      o << "Mean hadron squared pT vs. parton: " << mean_jt2 << std::endl;
      o << "Hadron pseudorapidity range: " << etamin << " to " << etamax << std::endl;
      o << "Minimum hadron pT (lab): " << min_hadron_pt << std::endl;
   }
   
   bool Configuration::Verify() {
      bool valid(true);
      // Check that the order is a valid range
      if(Order() < 0 or Order() > 1) {
         std::cerr << "Invalid order parameter:" <<
         "\n\tLO  = 0 (no pythia IS/FS interactions, LO FF)" <<
         "\n\tNLO = 1 (pythia IS/FS ON, NLO FF)" << std::endl;
         valid = false;
      } // if
      // Check that fragmentation functions that require a particular
      // order have that order.
      // Give warning but not error (valid still == true).
      if(Order() == LO and FragmentationFunction() == AKK) {
         std::cerr <<
         "WARNING: LO is requested but AKK FF's are only NLO. Calculation may not be self-consistent."
         << std::endl;
      } // if
      if(Order() == NLO and FragmentationFunction() == DSSC) {
         std::cout <<
         "WARNING: NLO is requested but Collins parametrization uses DSS LO. Order forced to 0 (LO)."
         << std::endl;
         order = LO;
      } // if
      return valid;
   }
   
} // namespace tppmc
